### MOVILIDAD EN ESTADOS DE FORMALIDAD Y COMERCIO INTERNACIONAL EN BOLIVIA ###
### Autor: Jos� Miguel Molina Fern�ndez ###

# Librer�as necesarias
library(haven)
library(labelled)
library(survey)
library(readxl)
library(car)
library(ggplot2)
library(extrafont)
library(margins)

# Importamos las bases de datos
# a) Encuestas Continuas de Empleo (se pueden descargar de / can be downloaded from: https://www.ine.gob.bo/index.php/banco/base-de-datos-sociales)
ece4_15 <- read_sav("D:/Mickey/Encuestas de Empleo/ECE_4T2015.sav")
for (i in 6:7) {
  for (j in 1:4) {
    ubica <- paste0("D:/Mickey/Encuestas de Empleo/ECE_", j, "T201", i, ".sav")
    name <- paste0("ece", j, "_1", i)
    assign(name, read_sav(ubica))
  }
}
rm(j) # limpiando
rm(name) # limpiando
rm(ubica) # limpiando

# b) Base de datos de crecimiento de exportaciones e importaciones del periodo
trade_growth <- read_excel("D:/Mickey/r4d trade labour markets informality/papers/determinantes de las transiciones entre formalidad e informalidad/data/trade_growth.xlsx", 
                           sheet = "Hoja3")

# Vectores para guardar cosas
att <- rep(NA, 6) # tasas de desgaste
personas1 <- rep(NA, 6) # personas en el panel en la primera entrevista

# Paneles (solo urbanos entran a paneles), tomando solamente primera y cuarta visitas (9 meses de distancia entre visitas)
panel4 <- rbind(ece4_15[ece4_15$panel==4,], ece3_16[ece3_16$panel==4,])
panel5 <- rbind(ece1_16[ece1_16$panel==5,], ece4_16[ece4_16$panel==5,])
panel6 <- rbind(ece2_16[ece2_16$panel==6,], ece1_17[ece1_17$panel==6,])
panel7 <- rbind(ece3_16[ece3_16$panel==7,], ece2_17[ece2_17$panel==7,])
panel8 <- rbind(ece4_16[ece4_16$panel==8,], ece3_17[ece3_17$panel==8,])
panel9 <- rbind(ece1_17[ece1_17$panel==9,], ece4_17[ece4_17$panel==9,])

# Nos quedamos solo con personas con empleo que:
# a) Tengan m�s de 14 a�os (edad m�nima de trabajo internacional)
# b) Est�n en el 99% mayor de la distribuci�n de ingresos reportada para las personas empleadas del rango de edad
# c) Trabajen en el sector privado
ldf <- list(panel4, panel5, panel6, panel7, panel8, panel9)
rm(panel4, panel5, panel6, panel7, panel8, panel9) # limpiando
names(ldf) <- paste0("panel", 4:9)
ldf <- lapply(ldf, function(df) {
 df <- df[df$condact==1 & !is.na(df$condact),]
 df <- df[df$s1_03a>=14,]
 df <- df[df$s2_22==3 | df$s2_22==4 | is.na(df$s2_22),]
 df$ylab[is.na(df$ylab)] <- -100 # previene problemas con NA's en siguiente subset
 df <- df[df$ylab==-100 | df$ylab>170.512,] # limite dado por menor percentil 1 observados en los paneles
 df$ylab[df$ylab==-100] <- NA
 return(df)
})

# Datos para calcular tasas de desgaste
personas1[1] <- nrow(ldf$panel4[ldf$panel4$trim==226,])
personas1[2] <- nrow(ldf$panel5[ldf$panel5$trim==227,])
personas1[3] <- nrow(ldf$panel6[ldf$panel6$trim==228,])
personas1[4] <- nrow(ldf$panel7[ldf$panel7$trim==229,])
personas1[5] <- nrow(ldf$panel8[ldf$panel8$trim==230,])
personas1[6] <- nrow(ldf$panel9[ldf$panel9$trim==231,])

# Indicadores antes de descartar a los que no aparecen en la cuarta entrevista
visita1 <- rbind(ldf$panel4[ldf$panel4$trim==226,], ldf$panel5[ldf$panel5$trim==227,], ldf$panel6[ldf$panel6$trim==228,], ldf$panel7[ldf$panel7$trim==229,], ldf$panel8[ldf$panel8$trim==230,], ldf$panel9[ldf$panel9$trim==231,])
sum(visita1$s1_02==1) / nrow(visita1)
sum(visita1$s1_10==1) / nrow(visita1)
table(visita1$s2_20) / nrow(visita1)
table(visita1$caeb_op) / nrow(visita1)
summary(visita1$yprilab)

# Nos quedamos solo con los que est�n en ambas visitas de la encuesta
ldf <- lapply(ldf, function(df) {
  duplis <- df$id_per_panel[duplicated(df$id_per_panel)==TRUE]
  df <- df[df$id_per_panel %in% duplis,] # mantiene solo a los que est�n en ambos periodos
  df <- copy_labels(ece4_15, df)
  return(df)
})

# Indicadores despues de descartar a los que no aparecen en la cuarta entrevista
visita1 <- rbind(ldf$panel4[ldf$panel4$trim==226,], ldf$panel5[ldf$panel5$trim==227,], ldf$panel6[ldf$panel6$trim==228,], ldf$panel7[ldf$panel7$trim==229,], ldf$panel8[ldf$panel8$trim==230,], ldf$panel9[ldf$panel9$trim==231,])
sum(visita1$s1_02==1) / nrow(visita1)
sum(visita1$s1_10==1) / nrow(visita1)
table(visita1$s2_20) / nrow(visita1)
table(visita1$caeb_op) / nrow(visita1)
summary(visita1$yprilab)

# Guardamos las tasas de desgaste
att[1] <- 1 - ((nrow(ldf$panel4)/2)/personas1[1])
att[2] <- 1 - ((nrow(ldf$panel5)/2)/personas1[2])
att[3] <- 1 - ((nrow(ldf$panel6)/2)/personas1[3])
att[4] <- 1 - ((nrow(ldf$panel7)/2)/personas1[4])
att[5] <- 1 - ((nrow(ldf$panel8)/2)/personas1[5])
att[6] <- 1 - ((nrow(ldf$panel9)/2)/personas1[6])

# Creamos variable de informalidad
ldf <- lapply(ldf, function(df) {
  df$informal <- ifelse((df$s2_20==1 | df$s2_20==2) & (df$s2_23==3 | (df$s2_23==4 & df$s2_64==2)) & df$s2_26<=5, 1,
                 ifelse(df$s2_20==3 & df$s2_23>=3, 1,
                 ifelse((df$s2_20==4 | df$s2_20==5) & df$s2_23>=3 & df$s2_26<=5, 1,
                 ifelse(df$s2_20>5, 1, 0))))
  df <- df[!is.na(df$informal), ] # gente para la que falta informacion (menos del 0.1% de cada panel)   
})     

# Volvemos a aplicar para omitir personas sin datos de informalidad en un periodo
ldf <- lapply(ldf, function(df) {
  duplis <- df$id_per_panel[duplicated(df$id_per_panel)==TRUE]
  df <- df[df$id_per_panel %in% duplis,] # mantiene solo a los que est�n en ambos periodos
  df <- copy_labels(ece4_15, df)
  return(df)
})

# Ultimos cambios:
# Trabajadores sin remuneraci�n de acuerdo a pregunta s2_20 tendr�n ingreso cero (para no perderlos en an�lisis de ingresos por ser NA)
# Cambiamos valor a extranjeros y personas que no se identifiquen como indigenas para que sean grupo de comparacion en regresiones
# Las personas que no respondieron seguro de salud no lo tienen, por tanto toman el valor correspondiente
# Creamos variable de educacion
ldf <- lapply(ldf, function(df) {
  df$yprilab[df$s2_20==7] <- 0
  df$s1_10[df$s1_10==2 | df$s1_10==3] <- 0
  df$s2_36a[is.na(df$s2_36a)] <- 2
  df$esup <- ifelse(df$s1_07a <=65, 0, 1)  # Tiene algun estudio superior?
  # la variable de a�os de escolaridad construida en la encuesta parece tener muchos problemas de NAs que no deber�an ser NAs, por eso se toma esta binaria
  return(df)
})

# Ordenamos los dataframes
ldf <- lapply(ldf, function(df) {
  df <- df[with(df, order(-gestion, -trimestre, id_per_panel)),]
})

# Guardamos los nuevos datasets, "_1" es el dataset de la primera visita y "_4" el de la segunda de cada panel
panel4_1 <- ldf$panel4[1:(nrow(ldf$panel4)/2),]
panel4_2 <- ldf$panel4[-(1:(nrow(ldf$panel4)/2)),]
panel5_1 <- ldf$panel5[1:(nrow(ldf$panel5)/2),]
panel5_2 <- ldf$panel5[-(1:(nrow(ldf$panel5)/2)),]
panel6_1 <- ldf$panel6[1:(nrow(ldf$panel6)/2),]
panel6_2 <- ldf$panel6[-(1:(nrow(ldf$panel6)/2)),]
panel7_1 <- ldf$panel7[1:(nrow(ldf$panel7)/2),]
panel7_2 <- ldf$panel7[-(1:(nrow(ldf$panel7)/2)),]
panel8_1 <- ldf$panel8[1:(nrow(ldf$panel8)/2),]
panel8_2 <- ldf$panel8[-(1:(nrow(ldf$panel8)/2)),]
panel9_1 <- ldf$panel9[1:(nrow(ldf$panel9)/2),]
panel9_2 <- ldf$panel9[-(1:(nrow(ldf$panel9)/2)),]

# Calculamos transiciones 1)formal-formal, 2) f-i, 3)i-f, 4)i-i
panel4_2$transicion <- ifelse(panel4_1$informal==0 & panel4_2$informal==0, 1,
                       ifelse(panel4_1$informal==0 & panel4_2$informal==1, 2,
                       ifelse(panel4_1$informal==1 & panel4_2$informal==0, 3, 4)))
panel5_2$transicion <- ifelse(panel5_1$informal==0 & panel5_2$informal==0, 1,
                       ifelse(panel5_1$informal==0 & panel5_2$informal==1, 2,
                       ifelse(panel5_1$informal==1 & panel5_2$informal==0, 3, 4)))
panel6_2$transicion <- ifelse(panel6_1$informal==0 & panel6_2$informal==0, 1,
                       ifelse(panel6_1$informal==0 & panel6_2$informal==1, 2,
                       ifelse(panel6_1$informal==1 & panel6_2$informal==0, 3, 4)))
panel7_2$transicion <- ifelse(panel7_1$informal==0 & panel7_2$informal==0, 1,
                       ifelse(panel7_1$informal==0 & panel7_2$informal==1, 2,
                       ifelse(panel7_1$informal==1 & panel7_2$informal==0, 3, 4)))
panel8_2$transicion <- ifelse(panel8_1$informal==0 & panel8_2$informal==0, 1,
                       ifelse(panel8_1$informal==0 & panel8_2$informal==1, 2,
                       ifelse(panel8_1$informal==1 & panel8_2$informal==0, 3, 4)))
panel9_2$transicion <- ifelse(panel9_1$informal==0 & panel9_2$informal==0, 1,
                       ifelse(panel9_1$informal==0 & panel9_2$informal==1, 2,
                       ifelse(panel9_1$informal==1 & panel9_2$informal==0, 3, 4)))

# Creamos una variable que nos diga si la persona ha cambiado de sector (ser� �til m�s tarde)
panel4_2$se_mueve <- ifelse(panel4_2$transicion==1 | panel4_2$transicion==4, 0, 1)
panel5_2$se_mueve <- ifelse(panel5_2$transicion==1 | panel5_2$transicion==4, 0, 1)
panel6_2$se_mueve <- ifelse(panel6_2$transicion==1 | panel6_2$transicion==4, 0, 1)
panel7_2$se_mueve <- ifelse(panel7_2$transicion==1 | panel7_2$transicion==4, 0, 1)
panel8_2$se_mueve <- ifelse(panel8_2$transicion==1 | panel8_2$transicion==4, 0, 1)
panel9_2$se_mueve <- ifelse(panel9_2$transicion==1 | panel9_2$transicion==4, 0, 1)

# Creamos dise�o muestral para cada panel
svypanel4 <- svydesign(ids=panel4_2$upm, strata=panel4_2$estrato, weights=panel4_2$fact_trim, data=panel4_2)
svypanel5 <- svydesign(ids=panel5_2$upm, strata=panel5_2$estrato, weights=panel5_2$fact_trim, data=panel5_2)
svypanel6 <- svydesign(ids=panel6_2$upm, strata=panel6_2$estrato, weights=panel6_2$fact_trim, data=panel6_2)
svypanel7 <- svydesign(ids=panel7_2$upm, strata=panel7_2$estrato, weights=panel7_2$fact_trim, data=panel7_2)
svypanel8 <- svydesign(ids=panel8_2$upm, strata=panel8_2$estrato, weights=panel8_2$fact_trim, data=panel8_2)
svypanel9 <- svydesign(ids=panel9_2$upm, strata=panel9_2$estrato, weights=panel9_2$fact_trim, data=panel9_2, nest=TRUE)
options(survey.lonely.psu="adjust") # correci�n de casos de un psu por estrato

# Calculamos el porcentaje de personas en cada estado de transici�n con intervalos de confianza
ldf <- list(svypanel4, svypanel5, svypanel6, svypanel7, svypanel8, svypanel9)
names(ldf) <- paste0("svypanel", 4:9)
lapply(ldf, function(df) {
  print("Transiciones del periodo")
  print(svyciprop(~I(transicion==1), df))
  print(svyciprop(~I(transicion==2), df))
  print(svyciprop(~I(transicion==3), df))
  print(svyciprop(~I(transicion==4), df))
  return(print(""))
})

# Encontramos la distribucion a la que converge todo
P = t(matrix(c(c(0.2082/(0.2082+0.0861), 0.0861/(0.2082+0.0861)),
               c(0.0841/(0.0841+0.6216), 0.6216/(0.0841+0.6216))),
             nrow=2)) # matriz de transici�n con medias geometricas de cada panel
A_basic <- t(diag(rep(1,nrow(P)))-P)
b_basic <- rep(0,nrow(P))
A_constr <- rbind(A_basic,rep(1,nrow(P)))
b_constr <- c(b_basic,1)
estacionaria <- t(solve(t(A_constr)%*%A_constr,t(A_constr)%*%b_constr))
estacionaria%*%P # Muestra la distribucion a la que se converge
rm(A_basic, b_basic, A_constr, b_constr) # limpiando

# Bases de datos de primera y cuarta visita del panel
visita1 <- rbind(panel4_1, panel5_1, panel6_1, panel7_1, panel8_1, panel9_1)
visita4 <- rbind(panel4_2, panel5_2, panel6_2, panel7_2, panel8_2, panel9_2)
visita1$transicion <- visita4$transicion
visita1$se_mueve <- visita4$se_mueve

# Cambios (o no) de afiiacion industrial de personas que cambian de estatus de formalidad
aux <- 0:20
print("Formal-Informal por cambio o no de actividad")
for (i in 1:21) {
  print(paste(sum((visita4$caeb_op[visita1$caeb_op!=visita4$caeb_op & visita4$transicion==2]==aux[i])*visita4$fact_trim[visita1$caeb_op!=visita4$caeb_op & visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2]), "actividad", i-1, sep="     "))
}
sum((visita4$caeb_op[visita1$caeb_op!=visita4$caeb_op & visita4$transicion==2]==99)*visita4$fact_trim[visita1$caeb_op!=visita4$caeb_op & visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
sum((visita1$caeb_op[visita4$transicion==2]==visita4$caeb_op[visita4$transicion==2]) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
print("Informal-formal por cambio o no de actividad")
for (i in 1:21) {
  print(paste(sum((visita4$caeb_op[visita1$caeb_op!=visita4$caeb_op & visita4$transicion==3]==aux[i])*visita4$fact_trim[visita1$caeb_op!=visita4$caeb_op & visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3]), "actividad", i-1, sep="     "))
}
sum((visita4$caeb_op[visita1$caeb_op!=visita4$caeb_op & visita4$transicion==3]==99)*visita4$fact_trim[visita1$caeb_op!=visita4$caeb_op & visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])
sum((visita1$caeb_op[visita4$transicion==3]==visita4$caeb_op[visita4$transicion==3]) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])
rm(aux, i) # limpiando

# Genero de personas que cambian de estatus de formalidad
print("Formal-Informal por genero")
sum((visita4$s1_02[visita4$transicion==2]==1) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
sum((visita4$s1_02[visita4$transicion==2]==2) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
print("Informal-formal por genero")
sum((visita4$s1_02[visita4$transicion==3]==1) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])
sum((visita4$s1_02[visita4$transicion==3]==2) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])

# Edad de personas que cambian de estatus de formalidad
print("Formal-informal por edad")
sum((visita4$s1_03a[visita4$transicion==2]<24) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
sum((visita4$s1_03a[visita4$transicion==2]>=24 & visita4$s1_03a[visita4$transicion==2]<34) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
sum((visita4$s1_03a[visita4$transicion==2]>=34 & visita4$s1_03a[visita4$transicion==2]<44) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
sum((visita4$s1_03a[visita4$transicion==2]>=44 & visita4$s1_03a[visita4$transicion==2]<54) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
sum((visita4$s1_03a[visita4$transicion==2]>=54 & visita4$s1_03a[visita4$transicion==2]<64) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
sum((visita4$s1_03a[visita4$transicion==2]>=64) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
print("Informal-formal por edad")
sum((visita4$s1_03a[visita4$transicion==3]<24) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])
sum((visita4$s1_03a[visita4$transicion==3]>=24 & visita4$s1_03a[visita4$transicion==3]<34) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])
sum((visita4$s1_03a[visita4$transicion==3]>=34 & visita4$s1_03a[visita4$transicion==3]<44) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])
sum((visita4$s1_03a[visita4$transicion==3]>=44 & visita4$s1_03a[visita4$transicion==3]<54) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])
sum((visita4$s1_03a[visita4$transicion==3]>=54 & visita4$s1_03a[visita4$transicion==3]<64) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])
sum((visita4$s1_03a[visita4$transicion==3]>=64) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])

# Identificaci�n ind�gena de personas que cambian de estatus de formalidad
print("Formal-Informal por identificaci�n ind�gena")
sum((visita4$s1_10[visita4$transicion==2]==1) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
print("Informal-formal por identificaci�n ind�gena")
sum((visita4$s1_10[visita4$transicion==3]==1) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])

# Personas con educacion superior que cambian de estatus de formalidad
print("Formal-Informal por nivel de educacion")
sum((visita4$esup[visita4$transicion==2]==1) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
print("Informal-formal por nivel de educacion")
sum((visita4$esup[visita4$transicion==3]==1) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])

# Promedio de ingresos laborales de personas que cambian de estatus de formalidad (antes)
print("Formal-informal por ingresos (antes)")
sum(visita1$yprilab[visita1$transicion==2] * visita1$fact_trim[visita1$transicion==2], na.rm=TRUE) / sum(visita1$fact_trim[visita1$transicion==2])
print("Informal-formal por ingresos (antes)")
sum(visita1$yprilab[visita1$transicion==3] * visita1$fact_trim[visita1$transicion==3], na.rm=TRUE) / sum(visita1$fact_trim[visita1$transicion==3])

# Promedio de ingresos laborales de personas que cambian de estatus de formalidad (despues)
print("Formal-informal por ingresos (despues)")
sum(visita4$yprilab[visita4$transicion==2] * visita4$fact_trim[visita4$transicion==2], na.rm=TRUE) / sum(visita4$fact_trim[visita4$transicion==2])
print("Informal-formal por ingresos (despues)")
sum(visita4$yprilab[visita4$transicion==3] * visita4$fact_trim[visita4$transicion==3], na.rm=TRUE) / sum(visita4$fact_trim[visita4$transicion==3])

# Promedio de horas trabajadas de personas que cambian de estatus de formalidad (antes)
print("Formal-informal por horas trabajadas (antes)")
sum(visita1$tothrs[visita1$transicion==2] * visita1$fact_trim[visita1$transicion==2], na.rm=TRUE) / sum(visita1$fact_trim[visita1$transicion==2])
print("Informal-formal por horas trabajadas (antes)")
sum(visita1$tothrs[visita1$transicion==3] * visita1$fact_trim[visita1$transicion==3], na.rm=TRUE) / sum(visita1$fact_trim[visita1$transicion==3])

# Promedio de horas trabajadas de personas que cambian de estatus de formalidad (despues)
print("Formal-informal por horas trabajadas (despues)")
sum(visita4$tothrs[visita4$transicion==2] * visita4$fact_trim[visita4$transicion==2], na.rm=TRUE) / sum(visita4$fact_trim[visita4$transicion==2])
print("Informal-formal por horas trabajadas (despues)")
sum(visita4$tothrs[visita4$transicion==3] * visita4$fact_trim[visita4$transicion==3], na.rm=TRUE) / sum(visita4$fact_trim[visita4$transicion==3])

# Cobertura de salud de personas que cambian de estatus de formalidad
print("Formal-informal por cobertura de salud (antes)")
sum((visita1$s2_36a[visita1$transicion==2]==1) * visita1$fact_trim[visita1$transicion==2]) / sum(visita1$fact_trim[visita1$transicion==2])
sum((visita1$s2_36a[visita1$transicion==2]==2) * visita1$fact_trim[visita1$transicion==2]) / sum(visita1$fact_trim[visita1$transicion==2])
print("Informal-formal por cobertura de salud (antes)")
sum((visita1$s2_36a[visita1$transicion==3]==1) * visita1$fact_trim[visita1$transicion==3]) / sum(visita1$fact_trim[visita1$transicion==3])
sum((visita1$s2_36a[visita1$transicion==3]==2) * visita1$fact_trim[visita1$transicion==3]) / sum(visita1$fact_trim[visita1$transicion==3])

# Cobertura de salud de personas que cambian de estatus de formalidad
print("Formal-informal por cobertura de salud (despues)")
sum((visita4$s2_36a[visita4$transicion==2]==1) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
sum((visita4$s2_36a[visita4$transicion==2]==2) * visita4$fact_trim[visita4$transicion==2]) / sum(visita4$fact_trim[visita4$transicion==2])
print("Informal-formal por cobertura de salud (despues)")
sum((visita4$s2_36a[visita4$transicion==3]==1) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])
sum((visita4$s2_36a[visita4$transicion==3]==2) * visita4$fact_trim[visita4$transicion==3]) / sum(visita4$fact_trim[visita4$transicion==3])

# Creamos variable de c�digo de CAEB a dos d�gitos transformado a CIIU rev. 3
visita1$s2_16acod <- ifelse(visita1$s2_16acod=="", NA, visita1$s2_16acod) # Las personas que no trabajaban estaban con un string vac�o en vez de valor de missing value (NA) en la base, lo cual hace dificil la identificacion de missing values
visita1$ciiu2 <- substr(visita1$s2_16acod, 1, 2)
indx <- match(visita1$ciiu2, c("01","02","03","05","06","07","08","10","11",as.character(12:25),"28","26","27","29","30","31","32","35","59"), nomatch=0)
visita1$ciiu2[indx!=0] <- c("01","02","05","10","11","13","14","15","15",as.character(16:23),"24","24",as.character(25:29),"30_32_33","31","34","35","36","36","40","92")[indx]
visita1$ciiu2[indx==0] <- "no_transables"
visita1$ciiu_trim <- paste0(visita1$ciiu2,"_", visita1$trim)
rm(indx) #limpiando

# Unimos la base de datos de la visita 1 con el crecimiento de importaciones desde a�o previo a visita 1 en actividad economica
visita1 <- merge(visita1, trade_growth, by="ciiu_trim", all.x=TRUE)

# Ajustes previos a regresiones de probabilidad de transicion 
visita1$upm_ae <- paste0(visita1$upm, "_", visita1$ciiu2) # para clusters
visita1 <- visita1[!is.na(visita1$yprilab),] # unica variable que sigue con NA's
visita1$yprilab <- visita1$yprilab / 100 # para mejor lectura de resultados

# Regresiones de probabilidad de transicion 
svyvis <- svydesign(ids=~upm_ae, strata=~estrato, weights=~fact_trim, data=visita1, nest=TRUE) # dise�o de la encuesta
svyvis1 <- subset(svyvis, transicion==1 | transicion==2) # mantiene informacion de clusters y estratos
svyvis2 <- subset(svyvis, transicion==3 | transicion==4) # mantiene informacion de clusters y estratos
print("Modelo de probabilidad de moverse del sector formal al informal")
model1 <- svyglm(as.factor(se_mueve) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(esup) + as.factor(s2_20) + as.factor(s2_36a) + as.factor(ciiu2) + as.factor(trim), family=binomial(link=logit), design=svyvis1)
summary(margins(model1, design=svyvis1)) # calcula efectos marginales promedio
print("Modelo de probabilidad de moverse del sector informal al formal")
model2 <- svyglm(as.factor(se_mueve) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(esup) + as.factor(s2_20) + as.factor(s2_36a) + as.factor(ciiu2) + as.factor(trim), family=binomial(link=logit), design=svyvis2)
summary(margins(model2, design=svyvis2)) # calcula efectos marginales promedio

# Significancia de variables de control con m�s de dos categor�as
model11 <- svyglm(as.factor(se_mueve) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(esup) + as.factor(s2_36a) + as.factor(ciiu2) + as.factor(trim), family=binomial(link=logit), design=svyvis1)
anova(model1, model11) 
model12 <- svyglm(as.factor(se_mueve) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(esup) + as.factor(s2_20) + as.factor(s2_36a) + as.factor(trim), family=binomial(link=logit), design=svyvis1)
anova(model1, model12)
model13 <- svyglm(as.factor(se_mueve) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(esup) + as.factor(s2_20) + as.factor(s2_36a) + as.factor(ciiu2), family=binomial(link=logit), design=svyvis1)
anova(model1, model13)
model21 <- svyglm(as.factor(se_mueve) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(esup) + as.factor(s2_36a) + as.factor(ciiu2) + as.factor(trim), family=binomial(link=logit), design=svyvis2)
anova(model2, model21)
model22 <- svyglm(as.factor(se_mueve) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(esup) + as.factor(s2_20) + as.factor(s2_36a) + as.factor(trim), family=binomial(link=logit), design=svyvis2)
anova(model2, model22)
model23 <- svyglm(as.factor(se_mueve) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(esup) + as.factor(s2_20) + as.factor(s2_36a) + as.factor(ciiu2), family=binomial(link=logit), design=svyvis2)
anova(model2, model23)

# Algunos diagnosticos de forma funcional y multicolinealidad
## probas1 <- predict(model1, type="response")
## probas2 <- predict(model2, type="response")
## logitp1 <- log(probas1/(1-probas1))
## logitp2 <- log(probas2/(1-probas2))
## plot(visita1$import_g[!is.na(visita1$tothrs) & !is.na(visita1$yprilab) & (visita1$transicion==1 | visita1$transicion==2)], logitp1)
## plot(visita1$import_g[!is.na(visita1$tothrs) & !is.na(visita1$yprilab) & (visita1$transicion==3 | visita1$transicion==4)], logitp2)
## plot(visita1$export_g[!is.na(visita1$tothrs) & !is.na(visita1$yprilab) & (visita1$transicion==1 | visita1$transicion==2)], logitp1)
## plot(visita1$export_g[!is.na(visita1$tothrs) & !is.na(visita1$yprilab) & (visita1$transicion==3 | visita1$transicion==4)], logitp2)
## plot(visita1$s1_03a[!is.na(visita1$tothrs) & !is.na(visita1$yprilab) & (visita1$transicion==1 | visita1$transicion==2)], logitp1)
## plot(visita1$s1_03a[!is.na(visita1$tothrs) & !is.na(visita1$yprilab) & (visita1$transicion==3 | visita1$transicion==4)], logitp2)
## plot(visita1$yprilab[!is.na(visita1$tothrs) & !is.na(visita1$yprilab) & (visita1$transicion==1 | visita1$transicion==2)], logitp1)
## plot(visita1$yprilab[!is.na(visita1$tothrs) & !is.na(visita1$yprilab) & (visita1$transicion==3 | visita1$transicion==4)], logitp2)
## plot(visita1$tothrs[!is.na(visita1$tothrs) & !is.na(visita1$yprilab) & (visita1$transicion==1 | visita1$transicion==2)], logitp1)
## plot(visita1$tothrs[!is.na(visita1$tothrs) & !is.na(visita1$yprilab) & (visita1$transicion==3 | visita1$transicion==4)], logitp2)
vif(model1) # multicolinealidad
vif(model2) # multicolinealidad

# Probabilidad de ser informal (Tabla 6) y signifancia de categ�ricas con m�s de dos categor�as
model3 <- svyglm(as.factor(informal) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(s2_36a) + as.factor(esup) + as.factor(s2_20) + as.factor(ciiu2) + as.factor(trim), family=binomial(link=logit), design=svyvis)
model4 <- svyglm(as.factor(informal) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(esup) + as.factor(s2_36a) + as.factor(s2_20) + as.factor(trim), family=binomial(link=logit), design=svyvis)
anova(model3, model4) 
model5 <- svyglm(as.factor(informal) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(esup) + as.factor(s2_36a) + as.factor(s2_20) +as.factor(ciiu2), family=binomial(link=logit), design=svyvis)
anova(model3, model5)
model6 <- svyglm(as.factor(informal) ~ import_g + export_g + as.factor(s1_02) + log(s1_03a) + as.factor(s1_10) + yprilab + tothrs + as.factor(esup) + as.factor(s2_36a) + as.factor(ciiu2) + as.factor(trim), family=binomial(link=logit), design=svyvis)
anova(model3, model6)
summary(margins(model3, design=svyvis)) # calculo de efectos marginales promedio
## stargazer(model3, ci=TRUE, report="vc*s", star.cutoffs = c(0.05, 0.01, 0.005)) # creacion de tabla base en formato LaTex

# Grafico 1 
distran2 <- data.frame(
  group = c("Trade and automobile repair", "Transportation and warehousing", "Construction", "Other industries", "Did not change industry affiliation"),
  value = c(8.77, 6.52, 5.19, 16.29, 63.23)
) # data

distran2 <- distran2[with(distran2, order(-value)),]
distran2$group <- factor(distran2$group, levels=distran2$group[c(3,4,5,2,1)])

graph1 <- ggplot(distran2, aes(x="", y=value, fill=group))+
  geom_bar(width = 1, stat = "identity") +
  scale_fill_brewer(palette="Set1") +
  theme_void() +
  theme(text = element_text(family="Garamond"),
        legend.title = element_blank(),
        legend.key.size = unit(0.3, "in"),
        legend.text = element_text(size=12)) +
  geom_text(aes(label = paste0(value, "%"),
                family = "Garamond",
                size = 10),
            position = position_stack(vjust = 0.5))

ggsave("D:\\Mickey\\r4d trade labour markets informality\\papers\\determinantes de las transiciones entre formalidad e informalidad\\graphs\\trans_graph1.png", graph1)

# Grafico 2 
distran3 <- data.frame(
  group = c("Manufacturing", "Trade and automobile repair", "Construction", "Other industries", "Did not change industry affiliation"),
  value = c(7.78, 7.30, 4.02, 17.07, 63.83)
)

distran3 <- distran3[with(distran3, order(-value)),]
distran3$group <- factor(distran3$group, levels=distran3$group[c(3,4,5,2,1)])

graph2 <- ggplot(distran3, aes(x="", y=value, fill=group))+
  geom_bar(width = 1, stat = "identity") +
  scale_fill_brewer(palette="Set3") +
  theme_void() +
  theme(text = element_text(family="Garamond"),
        legend.title = element_blank(),
        legend.key.size = unit(0.3, "in"),
        legend.text = element_text(size=12)) +
  geom_text(aes(label = paste0(format(round(value, 2), nsmall = 2), "%"),
                family = "Garamond",
                size = 10),
            position = position_stack(vjust = 0.5))

ggsave("D:\\Mickey\\r4d trade labour markets informality\\papers\\determinantes de las transiciones entre formalidad e informalidad\\graphs\\trans_graph2.png", graph2)

# Tabla 4
sum((visita1$informal==1 & visita1$ciiu2=="no_transables") * visita1$fact_trim) / sum(visita1$fact_trim[visita1$informal==1])
sum((visita1$informal==1 & visita1$ciiu2!="no_transables") * visita1$fact_trim) / sum(visita1$fact_trim[visita1$informal==1])
sum((visita1$informal==0 & visita1$ciiu2=="no_transables") * visita1$fact_trim) / sum(visita1$fact_trim[visita1$informal==0])
sum((visita1$informal==0 & visita1$ciiu2!="no_transables") * visita1$fact_trim) / sum(visita1$fact_trim[visita1$informal==0])

# Tabla 5
svy4_17 <- ece4_17
svy4_17$s2_16acod <- ifelse(svy4_17$s2_16acod=="", NA, svy4_17$s2_16acod) # Las personas que no trabajaban estaban con un string vac�o en vez de valor de missing value (NA) en la base, lo cual hace dificil la identificacion de missing values
svy4_17$ciiu2 <- substr(svy4_17$s2_16acod, 1, 2)
indx <- match(svy4_17$ciiu2, c("01","02","03","05","06","07","08","10","11",as.character(12:25),"28","26","27","29","30","31","32","35","59"), nomatch=0)
svy4_17$ciiu2[indx!=0] <- c("01","02","05","10","11","13","14","15","15",as.character(16:23),"24","24",as.character(25:29),"30_32_33","31","34","35","36","36","40","92")[indx]
svy4_17$ciiu2[indx==0] <- "no_transables"

svy4_17 <- svy4_17[svy4_17$condact==1 & !is.na(svy4_17$condact),]
svy4_17 <- svy4_17[svy4_17$s1_03a>=14,]
svy4_17 <- svy4_17[svy4_17$s2_22==3 | svy4_17$s2_22==4 | is.na(svy4_17$s2_22),]
svy4_17$ylab[is.na(svy4_17$ylab)] <- -100 # previene problemas con NA's en siguiente subset
svy4_17 <- svy4_17[svy4_17$ylab==-100 | svy4_17$ylab>170.512,] # limite dado por menor percentil 1 observados en los paneles
svy4_17$ylab[svy4_17$ylab==-100] <- NA
svy4_17$informal <- ifelse((svy4_17$s2_20==1 | svy4_17$s2_20==2) & (svy4_17$s2_23==3 | (svy4_17$s2_23==4 & svy4_17$s2_64==2)) & svy4_17$s2_26<=5, 1,
ifelse(svy4_17$s2_20==3 & svy4_17$s2_23>=3, 1,
ifelse((svy4_17$s2_20==4 | svy4_17$s2_20==5) & svy4_17$s2_23>=3 & svy4_17$s2_26<=5, 1,
ifelse(svy4_17$s2_20>5, 1, 0))))
svy4_17 <- svy4_17[!is.na(svy4_17$informal),]

sum((svy4_17$informal==1 & svy4_17$ciiu2=="no_transables") * svy4_17$fact_trim) / sum(svy4_17$fact_trim[svy4_17$informal==1])
sum((svy4_17$informal==1 & svy4_17$ciiu2!="no_transables") * svy4_17$fact_trim) / sum(svy4_17$fact_trim[svy4_17$informal==1])
sum((svy4_17$informal==0 & svy4_17$ciiu2=="no_transables") * svy4_17$fact_trim) / sum(svy4_17$fact_trim[svy4_17$informal==0])
sum((svy4_17$informal==0 & svy4_17$ciiu2!="no_transables") * svy4_17$fact_trim) / sum(svy4_17$fact_trim[svy4_17$informal==0])